knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(message = FALSE)
knitr::opts_chunk$set(warning = FALSE)
knitr::opts_chunk$set(cache = TRUE)
library(tidyverse)
library(cowplot)
library(broom)
library(readxl)
library(janitor)
library(plotrix)
library(here)
library(growthTools)
library(rootSolve)
library(nls.multstart)
library(nlstools)
Read in data
treatments <- read_excel(here("data-general", "ChlamEE_Treatments_JB.xlsx")) %>%
clean_names() %>%
mutate(treatment = ifelse(is.na(treatment), "none", treatment)) %>%
filter(population != "cc1629")
light_data <- read_csv(here("data-processed", "light-rstar-rfus-time.csv"))
ldata <- light_data %>%
mutate(percentage = ifelse(percentage == "0.5-0.7", "0.6", percentage)) %>%
mutate(percentage = as.numeric(percentage)) %>%
mutate(percentage = percentage/100) %>%
mutate(light = percentage*250)
Estimate growth rates, using piecewise regression (either exponential, lagged growth, saturated or lagged then saturated). Here we are choosing between moddels using AICc.
growth_rates_n_AICc <- ldata %>%
filter(population != "COMBO") %>%
mutate(ln.fluor = log(RFU)) %>%
group_by(well_plate) %>%
do(grs = get.growth.rate(x = .$days, y = .$ln.fluor,id = .$well_plate, plot.best.Q = F))
Pull out the things we want
growth_sum_n_AICc <- growth_rates_n_AICc %>%
summarise(well_plate, mu = grs$best.slope, n_obs = grs$best.model.slope.n,
slope_r2 = grs$best.model.slope.r2,
best.model_r2 = grs$best.model.rsqr,
best.model = grs$best.model, best.se = grs$best.se,
contents = grs$best.model.contents)
key <- ldata %>%
select(well_plate, population, light) %>%
distinct(well_plate, .keep_all = TRUE)
AICc_growth_rates <- growth_sum_n_AICc %>%
select(-contents) %>%
mutate(IC_method = "AICc")
AICc_growth_rates2 <- left_join(AICc_growth_rates, key, by = "well_plate")
# write_csv(AICc_growth_rates2, here("data-processed", "nitrate_exp_growth_w_growthtools_AICc.csv"))
exp_params <- growth_sum_n_AICc %>%
unnest(contents %>% map(tidy, .id = "number")) ## pull out the slopes and intercepts etc.
exp_params_aug <- growth_sum_n_AICc %>%
unnest(contents %>% map(augment, .id = "number")) %>% ### now add the fitted values
rename(days = x)
exp_wide <- exp_params %>%
spread(key = term, value = estimate)
all_preds <- left_join(exp_params_aug, exp_wide)
all_preds2 <- left_join(all_preds, key, by = "well_plate") %>%
group_by(well_plate) %>%
mutate(B1 = mean(B1, na.rm = TRUE)) %>%
mutate(B2 = mean(B2, na.rm = TRUE)) %>% ### here we do some wrangling to get the colour coding right for our plots
mutate(exponential = case_when(best.model == "gr" ~ "yes",
best.model == "gr.sat" & days < B1 ~ "yes",
best.model == "gr.lag" & days < B1 ~ "yes",
best.model == "gr.lagsat" & days < B2 & days > B1 ~ "yes",
TRUE ~ "no"))
Plot all the time series, here colour-coded by which model was the best fit, and the shape of the points corresponds to whether the point was in the exponential phase
all_preds2 %>%
ggplot(aes(x = days, y = .fitted, group = well_plate, color = best.model)) + geom_line() +
geom_point(aes(x = days, y = y, shape = exponential)) +
facet_grid(light ~ population) + ylab("Ln(RFU)") +xlab("Days")
ggsave(here("figures", "light_time_series_growth_tools.png"), width = 45, height = 20)
Fit Eilers-Peeters model
epcurve <- function(light, alpha, eopt, ps, m){
res <- (light/((1/(alpha *
eopt^2)) * light^2 + (1/ps - 2/(alpha * eopt)) * light + (1/alpha))) + h
res
}
AICc_growth_rates3 <- left_join(AICc_growth_rates2, treatments, by = "population")
ep_fits_i_m <- AICc_growth_rates3 %>%
rename(r_estimate = mu) %>%
distinct(r_estimate, light, population, .keep_all = TRUE) %>%
group_by(population, treatment, ancestor_id) %>%
nest() %>%
mutate(fit = purrr::map(data, ~ nls_multstart(r_estimate ~ (light/((1/(alpha *
eopt^2)) * light^2 + (1/ps - 2/(alpha * eopt)) * light + (1/alpha)))-0.1,
data = .x,
iter = 500,
start_lower = c(alpha = 0.2, eopt = 100, ps = 1),
start_upper = c(alpha = 0.4, eopt = 200, ps = 2),
supp_errors = 'N',
na.action = na.omit,
lower = c(alpha = 0, eopt = 80, ps = 0.1),
upper = c(alpha = 4, eopt = 400, ps = 5),
control = nls.control(maxiter=1000, minFactor=1/204800000))))
fits_many <- ep_fits_i_m
params <- fits_many %>%
filter(fit != "NULL") %>%
unnest(fit %>% map(tidy))
info <- fits_many %>%
unnest(fit %>% map(glance))
ep_wide <- params %>%
select(population, term, estimate) %>%
spread(key = term, value = estimate, 2:3)
CI <- fits_many %>%
filter(fit != "NULL") %>%
unnest(fit %>% map(~ confint2(.x) %>%
data.frame() %>%
dplyr::rename(., conf.low = X2.5.., conf.high = X97.5..))) %>%
group_by(., population) %>%
mutate(., term = c('alpha', 'eopt', 'ps')) %>%
ungroup()
# merge parameters and CI estimates
params <- merge(params, CI, by = intersect(names(params), names(CI)))
write_csv(params, here("data-processed", "EP-with-mortality-parameters.csv"))
# get predictions
preds <- fits_many %>%
unnest(fit %>% map(augment))
new_preds <- AICc_growth_rates3 %>%
rename(estimate = mu) %>%
distinct(estimate, light, population, .keep_all = TRUE) %>%
group_by(population, treatment, ancestor_id) %>%
mutate(r_estimate= estimate) %>%
do(., data.frame(light = seq(min(.$light), max(.$light), length.out = 150), stringsAsFactors = FALSE))
# max and min for each curve
p_growth3 <- AICc_growth_rates3 %>%
rename(estimate = mu) %>%
distinct(estimate, light, population, .keep_all = TRUE)
max_min <- dplyr::group_by(p_growth3, population, treatment, ancestor_id) %>%
summarise(., min_light = min(light), max_light = max(light)) %>%
ungroup()
# create new predictions
preds2 <- fits_many %>%
unnest(fit %>% map(augment, newdata = new_preds)) %>%
merge(., max_min) %>%
group_by(., population) %>%
filter(., light > unique(min_light) & light < unique(max_light)) %>%
dplyr::rename(., r_estimate = .fitted) %>%
ungroup()
Visualize the fits of the EP model to the growth rate data
ggplot() +
geom_point(aes(light, estimate), size = 2, data = p_growth3) +
geom_line(aes(light, r_estimate, group = population), data = preds2) +
facet_grid(ancestor_id ~ treatment, scales = "free") +
ylab('Growth rate (/day)') +
xlab('Irradiance')
Now find I* from root of EP, set mortality rate (respiratory losses) to 0.1 (as in m = 0.1 for the nitrate analysis)
epcurve_m<-function(light, alpha, eopt, ps){
res <- (light/((1/(alpha *
eopt^2)) * light^2 + (1/ps - 2/(alpha * eopt)) * light + (1/alpha)))-0.1
res
}
get_rstar_i <- function(df){
roots <- uniroot.all(function(x) epcurve_m(x, alpha = df$estimate[[1]], eopt = df$estimate[[2]],
ps = df$estimate[[3]]), interval = c(0, 10))
rstar <- roots[[1]]
rstar2 <- as.data.frame(rstar)
return(rstar)
}
ep_fits_split <- params %>%
split(.$population)
rstars_ep <- ep_fits_split %>%
map(get_rstar_i) %>%
unlist() %>%
as_data_frame() %>%
dplyr::rename(rstar = value) %>%
mutate(population = rownames(.))
rstars2_ep <- left_join(rstars_ep, treatments, by = "population")
Plot the resulting I*
rstars2_ep %>%
group_by(treatment) %>%
summarise_each(funs(mean, std.error), rstar) %>%
ggplot(aes(x = reorder(treatment, mean), y = mean)) + geom_point() +
geom_errorbar(aes(ymin = mean - std.error, ymax = mean + std.error), width = 0.1) +
# geom_boxplot(aes(x = treatment, y = rstar), color= "black", data = filter(rstars2_ep, rstar > 0)) +
geom_point(aes(x = reorder(treatment, rstar), y = rstar, color = ancestor_id), data =rstars2_ep, size = 2) + ylab("I* (umols/m2/s)") +
xlab("Selection treatment")